namespace Dune {
  namespace Compose {
    template<typename EQ, typename INIT, typename BOUN>
    void assigngrid(EQ& eq, INIT& init, BOUN& boun, Dune::ParameterTree& ptree) {
      const int dim = EQ::dim;
      const int nvar = EQ::nvar;
      const int refine = ptree.get<int>("grid.refinement");
      const int degree = ptree.get<int>("fem.degree");
      const std::string manager = ptree.get("grid.manager","yasp");
      /**
       * select grid
       */
      // 1D YASP for 1D finite volume (FV) method test
      if (manager=="yasp" && dim==1) {
        const int dim = 1; 
        typedef Dune::YaspGrid<1> Grid;
        typedef Grid::ctype DF;
        Dune::FieldVector<DF,1> L;
        L[0] = ptree.get("grid.structured.LX",(double)1.0);
        std::array<int,1> N;
        N[0] = ptree.get("grid.structured.NX",(int)16);
        std::bitset<1> B(false);
        int overlap=1;
        std::shared_ptr<Grid> gridp = std::shared_ptr<Grid>
          (new Grid(L,N,B,overlap,Dune::MPIHelper::getCollectiveCommunication()));
        gridp->refineOptions(false); // keep overlap in cells
        // gridp->globalRefine(refine);
        typedef Grid::LeafGridView GV;
        GV gv=gridp->leafGridView();
        using FLUX = Dune::Flux::FluxDeduce<EQ,GV::Grid::ctype>;
        FLUX flux(eq,ptree); 
        if (degree==0) {
          std::cout << "assgin grid: YASP, dim=2, deg=0" << std::endl;
          typedef Dune::PDELab::QkDGLocalFiniteElementMap
            <GV::Grid::ctype,double,0,1,Dune::PDELab::QkDGBasisPolynomial::legendre> FEM;
          FEM fem;
          Dune::Compose::driver(eq,flux,init,boun,gv,fem,ptree);
        }
        if (degree==1) {
          std::cout << "assgin grid: YASP, dim=2, deg=1" << std::endl;
          typedef Dune::PDELab::QkDGLocalFiniteElementMap
            <GV::Grid::ctype,double,1,1,Dune::PDELab::QkDGBasisPolynomial::legendre> FEM;
          FEM fem;
          Dune::Compose::driver(eq,flux,init,boun,gv,fem,ptree);
        }
        if (degree==2) {
          std::cout << "assgin grid: YASP, dim=2, deg=1" << std::endl;
          typedef Dune::PDELab::QkDGLocalFiniteElementMap
            <GV::Grid::ctype,double,2,1,Dune::PDELab::QkDGBasisPolynomial::legendre> FEM;
          FEM fem;
          Dune::Compose::driver(eq,flux,init,boun,gv,fem,ptree);
        }
      }

      // 2D YASP for 2D FV and Discontinous Galkin (DG) method
      if (manager=="yasp" && dim==2) {
        const int dim = 2;
        typedef Dune::YaspGrid<2> Grid;
        typedef Grid::ctype DF;
        Dune::FieldVector<DF,2> L;
        L[0] = ptree.get("grid.structured.LX",(double)1.0);
        L[1] = ptree.get("grid.structured.LY",(double)1.0);
        std::array<int,2> N;
        N[0] = ptree.get("grid.structured.NX",(int)16);
        N[1] = ptree.get("grid.structured.NY",(int)16);
        std::bitset<2> B(false);
        int overlap=1;
        std::shared_ptr<Grid> gridp = std::shared_ptr<Grid>
          (new Grid(L,N,B,overlap,Dune::MPIHelper::getCollectiveCommunication()));
        gridp->refineOptions(false); // keep overlap in cells
        // std::cout << "refine:" << refine << std::endl;
        gridp->globalRefine(refine);
        typedef Grid::LeafGridView GV;
        GV gv=gridp->leafGridView();
  
        using FLUX = Dune::Flux::FluxDeduce<EQ,GV::Grid::ctype>;
        FLUX flux(eq,ptree); 
  
        if (degree==0) {
          std::cout << "assgin grid: YASP, dim=2, deg=0" << std::endl;
          typedef Dune::PDELab::QkDGLocalFiniteElementMap
            <GV::Grid::ctype,double,0,2,Dune::PDELab::QkDGBasisPolynomial::legendre> FEM;
          FEM fem;
          Dune::Compose::driver(eq,flux,init,boun,gv,fem,ptree);
        }
        if (degree==2) {
          std::cout << "assgin grid: YASP, dim=2, deg=1" << std::endl;
          typedef Dune::PDELab::QkDGLocalFiniteElementMap
            <GV::Grid::ctype,double,2,2,Dune::PDELab::QkDGBasisPolynomial::legendre> FEM;
          FEM fem;
          Dune::Compose::driver(eq,flux,init,boun,gv,fem,ptree);
        }
      }
  
      // 4D YASP for 3D FV and DG method
      if (manager=="yasp" && dim==3) {
        const int dim = 3;
        typedef Dune::YaspGrid<3> Grid;
        typedef Grid::ctype DF;
        Dune::FieldVector<DF,3> L;
        L[0] = ptree.get("grid.structured.LX",(double)1.0);
        L[1] = ptree.get("grid.structured.LY",(double)1.0);
        L[2] = ptree.get("grid.structured.LZ",(double)1.0);
        std::array<int,3> N;
        N[0] = ptree.get("grid.structured.NX",(int)16);
        N[1] = ptree.get("grid.structured.NY",(int)16);
        N[2] = ptree.get("grid.structured.NZ",(int)16);
        std::bitset<3> B(false);
        int overlap=1;
        std::shared_ptr<Grid> gridp = std::shared_ptr<Grid>
          (new Grid(L,N,B,overlap,Dune::MPIHelper::getCollectiveCommunication()));
        gridp->refineOptions(false); // keep overlap in cells
        // gridp->globalRefine(refine);
        typedef Grid::LeafGridView GV;
        GV gv=gridp->leafGridView();

        using FLUX = Dune::Flux::FluxDeduce<EQ,GV::Grid::ctype>;
        FLUX flux(eq,ptree); 

        if (degree==0) {
          typedef Dune::PDELab::QkDGLocalFiniteElementMap
            <GV::Grid::ctype,double,0,3,Dune::PDELab::QkDGBasisPolynomial::legendre> FEM;
          FEM fem;
          Dune::Compose::driver(eq,flux,init,boun,gv,fem,ptree);
        }
        if (degree==1) {
          typedef Dune::PDELab::QkDGLocalFiniteElementMap
            <GV::Grid::ctype,double,1,3,Dune::PDELab::QkDGBasisPolynomial::legendre> FEM;
          FEM fem;
          // Dune::Compose::driver(eq,flux,init,boun,gv,fem,ptree);
        }        
      }

      // 2D UG for 2D FV and DG method
      if (manager=="ug" && dim==2) {
      }

      // 3D UG for 3D FV and DG method
      if (manager=="ug" && dim==3) {
      }






    } // driver
  }  //Compose
}  // Dune 
